arXiv:1504.02140vl [astro-ph.SR] 8 Apr 2015 


Accepted to ApJ: April 6, 2015 

Preprint typeset using style emulateapj v. 5/2/11 


PRESTELLAR CORE EORMATION, EVOLUTION, AND ACCRETION EROM GRAVITATIONAL 
FRAGMENTATION IN TURBULENT CONVERGING FLOWS 

Munan Gong and Eve C. Ostriker 

Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544, USA 

Accepted to ApJ: April 6, 2015 

ABSTRACT 

We investigate prestellar core formation and accretion based on three-dimensional hydrodynamic 
simulations. Our simulations represent local ~ Ipc regions within giant molecular clouds where a 
supersonic turbulent flow converges, triggering star formation in the post-shock layer. We include 
turbulence and self-gravity, applying sink particle techniques, and explore a range of inflow Mach 
number M. = 2 — 16. Two sets of cores are identified and compared: ti-cores are identified of a 
time snapshot in each simulation, representing dense structures in a single cloud map; tcoii-cores are 
identified at their individual time of collapse, representing the initial mass reservoir for accretion. We 
find that cores and filaments form and evolve at the same time. At the stage of core collapse, there is 
a well-defined, converged characteristic mass for isothermal fragmentation that is comparable to the 
critical Bonner-Ebert mass at the post-shock pressure. The core mass functions (CMFs) of Uoii-cores 
show a deficit of high-mass cores (> 7Mq) compared to the observed stellar initial mass function 
(IMF). However, the CMFs of ti-cores are similar to the observed CMFs and include many low-mass 
cores that are gravitationally stable. The difference between U-cores and tcoii-cores suggests that the 
full sample from observed CMFs may not evolve into protostars. Individual sink particles accrete 
at a roughly constant rate throughout the simulations, gaining one tcoircore mass per free-fall time 
even after the initial mass reservoir is accreted. High-mass sinks gain proportionally more mass at 
late times than low-mass sinks. There are outbursts in accretion rates, resulting from clumpy density 
structures falling into the sinks. 

Subject headings: ISM: structure - methods: numerical - stars: formation — stars: luminosity func¬ 
tion, mass function - turbulence 


1. INTRODUCTION 

Dense molecular cores form in the last step of giant 
molecular cloud (GMC) fragmentation and a re the cra- 


and are often associated with dense filaments seen 
in dust continuum and molecular lines. Observations 
show many pressure-confined gravitationally stable cores 


Ostriker 

2007 

Andre et al. 

2014 


__ __ These cores provide 

the initial mass reservoir and set the local environment 
of star formation. Understanding core formation and the 
core mass function (CMF) is important for characterizing 
the initial conditions of star and planet formation, and 
explaining the origi n of the stellar ini tial mass function 
(IMF, see review of Offner et al.||2014 ). 


gravitationallv unstable ( 

Lada et al. 

2008 

Andre et al. 

2010 

Schisano et al.|2014 

). High resolution observations 


ally bound cores are m ainly embedded along thermally 


super-critical filaments (lAndre et al. 12010 

Molinari et al. 

2010 Schisano et al. 

2014D, i.e., t 

rose with mass per 

unit length exceeding 2c^lG (Ostriker 

1964). This is 


Cores emerge in dense regions of GMCs, the dynamics 
of which ar e governed by tur bulence, gravity, and mag¬ 
netic fields (Dobbs et al.f2014 1. Although supersonic tur- 


that self-gravitatin g filaments are unstable to longitudi- 


bulence provides support tor GMCs on global scales, it ^992 1997 1, which would allow them to fragrnent into 


nal perturb ations ( Nagasawa| 1987 Inutsuka & Miyama 


gravitational collapse can rapidly occur 

(jBallesteros- 

Paredes et al. 2007 

McKee & Ostriker 2007 

). Strong tur- 


and aid in the creation of magnetically supercritica l cores 


1 Nakamura & Li|| 

2005 

2008 Kudoh & Basu||2008 

2011 

Chen & Ostriker| 

2012 

. However, in the situation where 

the large-scale Gi 

VIC is magnetically supercritical, simu- 


lations show that anisotropic core formation via magnetic 
field-aligned flows results in similar core m asses to the un- 


2014). In this paper. 


magnetized case (Chen & Ostriker 

we focus on core formation and subsequent collapse and 
accretion using turbulent, unmagnetized simulations. 
Cores are at the bottom of the hierarchy in a GMC, 

munan@princeton.edu, eco@astro.princeton.edu 


dense cores. However, the most unstable scale is several 
times the filament diameter, which would produce cores 
more massive than observed. Also, the two-step scenario 
in which filaments first form and then subsequently frag¬ 
ment into cores can be altered with the presence of large- 
amplitude turbulence, because the nonlinear perturba¬ 
tions produced by turbulence allow structures at many 
scales to grow simultaneously. With numerical simula¬ 
tions, it is possible to follow the temporal evolution of 
developing structures, and simulations with turbulence 
show that filaments and cores develop at the same time 
rather than in sequence ( Gong fc Ostriker|20Il Van Loo| 


let al.||2014||Chen fe Ustriker||2L)14[ see also p.i.l| in this 

paperj^ 

Observationally, cores are divided into two categories: 
starless cores, and protostellar cores/Class 0 objects (HI 
















































































2 


Francesco et al. 2007 Andre et al. 2009). Among the 
starless cores, a lurther classihcation is otprestellar cores, 
which are strongly gravitationally bound and centrally 
concentrated; the balance of other starless cores are the 
opposite: they appear to be largely pressure-confined 
rather than strongly gravitationally bound, and are not 
strongly centrally concentrated. Protostellar cores/Class 
0 objects are gravitationally bound and contain em¬ 
bedded protostars/young stellar objects (YSOs). It is 
generally thought that these categories represent cores 
at different evolutionary stages: pressure-confined star¬ 
less cores gain mass and turn into gravitationally bound 
prestellar cores, which finally collapse and become pro¬ 
tostellar cores. However, in practice, prestellar cores are 
often not distinguished from other starless cores in ob¬ 
servations. Some starless cores may even be transient 
structures that will never collapse to form protostars. 
Numerical simulation made it possible to check whether 
observed pressure-confined cores are likely to become un¬ 
stable at a later time. 

Many theoretical studies have investigated individual 
co re evolution, b eginning with t he se mi-analytic al work 


teristic timescale for the core collapse and the envelope 
infall stages is the free-fall time 






- 1/2 


• ( 3 ) 


This is similar to the lifet ime of prestellar an d protostel¬ 
lar cores in observations ( Andre et al.||2014 ). 

Empirical measurement of ClVlfi's from extinction and 
continuum mapping suggest the mass distribution has 
a similar shape to the IMF, but with the charact eristic 

core mass shifted upward by a factor of _ ^_ 

et al. 1 12007 1 |Nutter & Ward-Thompson||2007| |Konyves 
et al.|2010p . However, observed CMFs often include tran- 


sient structures that may not collapse later, and the core 
masses obtained from observations can be dependent on 
the spec ific core finding algorithm applied (e. g., gauss- 


clump s (Stutzki fc Guesten 19901); clumpGnd (Williams 


et al.||19gl ); getsources ( [Meri'shcnikov et al.|[20^ )). Al- 


of Larson (1969), Pension (1969), and Shu (1977). Based 
on isothermal simulations wit 


flows, Gong & Ostriker 
core evolution: core bui. 


1 supersonic converging 2009 Gong & Ostriker 


though core identification based directly on surface den¬ 
sity maps is somewhat ambiguous due to the arbitrary 
definition of the core boundary, use of the gravitational 
potential removes some of th is ambiguity (Smith et al. 


(2009) identified four stages of 
-ding, core collapse, envelope in- 


2011 )^_Analytical theories of the 


GIVIF (see review of Uftner et al.||2014 ) focus on gravita- 


fall, and late accretion. In the core building stage, the 
core gains mass from its environment, gradually becom¬ 
ing more centrally concentrated. For the spherical case, 
shock-bounded cores collapse when the central-to-edge 
pressure ratio is si milar to that of a critical Bo nnor-Ebert 
equilibrium sphere (Ebert 1955 Bonnor 1956), which has 


mass 


Mbe = 1-2 


(G3Pedge)l/2 
= 2.3Mo ( 

° V 10^cm“3/ 


= 1.9 


(G3p)l/2 
— 1/2 / J' \ 3/2 

[imj 


( 1 ) 


Here Hedge is external pressure at the edge of the core; 
p is the mean core density; T the core’s temperature; 
and Cs = {kT/pLY^"^ is the internal sound speed in the 
core. After the core collapses, the internal density struc- 
ture approaches the well-kno wn Larson-Penston solution 
( Larson|[l969 Penston||1969 hereafter LP) 


tionally bound cores. In this work, we identify two sets 
of cores: The first set is comparable to structures seen in 
a GMC at a given instant in time, and includes a range 
of evolutionary stages (both weakly and strongly bound). 
The second set represents cores seen at a fixed evolution¬ 
ary state, the time of protostar formation. We address 
the differences between the two sets, and implications for 
interpretation of observations. 

Even if every observed core were to collapse, evolution 
after core formation could further change the relation be¬ 
tween the GMF and IMF. At a minimum, fragmentation 
into binaries and multiples due to rapid rotation, and 
mass loss due to protostellar outflows, complicates any 
mapping from the GMF to the IMF. In addition, a con¬ 
cern in the theoretical literature has been that there may 
be no preferred scale for isothermal fragmentation such 
that gravitational fragmentation in GMGs must depend 


crucially on non-isothermal effects (Martel et al. 2006 


Krumholz 

2014 

). 'Lhis concern arises in part because an 

isotherma 

filament with mass per unit length greater 


than 2Cs/G can collapse towards it s axis without frag- 


PLp(?') = 


SMci 
47rGr2 ’ 


( 2 ) 


menting, leading to infinite density (Inutsuka & Miyama 


1997), and in part because homologous collapse of a 


where r is the distance from the core center. This asymp¬ 
totic state is reached regardless of the dynamical evo¬ 
lution le ading up to collapse (see disc ussion and refer¬ 
ence s in |Goi^ fc Ostriker] |20()9[ |201H hereafter |GO09| 
and GOlip ; in particular, cores formed by supersonic 
flows approach the LP solution during collapse. Collapse 
leads to creation of an opaque core and then a p roto¬ 


sphere would lead to ever smaller Jeans masses and pos¬ 
sible su b-fragmenta tion until stopped by non-isothermal 
effects ( Hoyle|[l9^ ). 


In fact, core collapse is non-homologous, such t h at the 


star at the center, with a surrounding rotating disk (Shu 
et al.||19^). Subsequently, the dense envelope falls onto 


the protostar via an inside-out gravitati o nally-in duced 
rarefaction wave (|Shu||1977| |Hunter||1977l) . |GO09| define 
the late accretion stage as the period after the rarefac- 
tion wave has propagated beyond the region where the 
LP profile holds, and most of the original core mass has 
fallen into the center (the protostar/disk). The charac¬ 


hierarchical fragmentation envisioned by Hoyle (1953) 
does not occur. However, disks that form around proto- 
stars are known to undergo excessive fragmentation and 
overproduce brown d warfs if stel lar heating is not in¬ 
cluded in simulations (Bate 2009). Here, we argue that 
the problem of excessive disk fragmentation in isother¬ 
mal models can be distinguished form whether there is a 
well-defined GMF from gravoturbulent fragmentation of 
an isothermal cloud. We shall show that the character¬ 
istic core mass converges with increasing resolution (and 
threshold of sink particle creation) in our simulations. 
In reality, of course, GMGs and cores are not isothermal. 
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but have temperature variations between ~ 5 —lOK, with 
observationally inferred te mperatures dropping toward 
the center of the cores (e.g., Marsh et al.|2014 Roy et al. 
20141. However, the isothermal approximatiori is an ad- 
equate first limit from cloud to core scales. Radiation 
feedback and significant non-isothermal effects become 
important subsequent to core collapse, in the accretion 
disks around protostars. 

In this paper, we present our results on prestellar 
core formation and collapse based on a large suite of 
three-dimensional numerical simulations. Our simula¬ 
tions model a local (~ Ipc) region of GMCs where a 
supersonic turbulent flow converges, inducing star forma¬ 
tion in the post-shock layer. We adopt the same model 
of a supersonic and isother mal con verging flow w ith tur¬ 
bulent perturbations as in |G011[ As in |G011[ we see 
filamentary structures form at the same time as dense 
cores. While the models of IGO 111 were halted when the 
first core collap sed, here we apply the sink particle im¬ 
plementation of Gong & Ostriker (2013) to trace the evo¬ 
lution at later times. We also explore a wider range of 
inflow Mach numbers (At = 2 — 16). Using simulations 
with varying numerical resolution, we show that the char¬ 
acteristic core mass at the time of collapse is well defined. 
We measure the statistical properties of cores, especially 
the CMF, and compare to the observed GMF and IMF. 
We are also able to trace the accretion rate of individual 
sink particles, and explore its time variation and depen¬ 
dence on core masses. 

The outline of this paper is as follows. We describe our 
numerical methods and model parameters in 1 ^ In 
we present our results of overall evo lutio n (q3T I, conver¬ 
gence of isot herm al fragmentation (( |3.2[ ), core prop erties 
and CMF ((|3.3|), and sink particle accretion ( |3.4[ ). Fi¬ 
nally summarizes our conclusions. 

2. METHODS AND PARAMETERS 

2.1. Numerical Methods 

2.1.1. Basic Equations and Algorithms 

The numerical simulations pre sented in this paper are 
conducted with th e Athena code (|Stone et al. 2008| Stone, 
fc Gardiner|2009 1, on a three dimensional (31J) Cartesian 
grid, using the van Leer (MUSCL-Hancock type) integra¬ 
tor, the HLLC Riemann solver, and second-order spa¬ 
tial reconstruction. The code solves for the self-gravity 
in planar geometry, with open boundary condition in 
the vertical (z) direction, and periodic boundary con¬ 
ditions in the horizontal (i and y) directions, using the 
fast Fourier transformation (FFT) method developed by 
Koyama & Ostriker (|2009|). Sink particles as imple- 
Ustrike] 


mented by Gong & Ostriker are included, and 


allow us to track the system evolution after individual 
cores collapse. 

The equations solved for the gas are the conservation 
of mass and momentum: 


§( + V.(pv)=0, 

5v VP 

—-PvVv =-V(d>-P$sp), 

at p 

and the Poisson equation: 

= inGp, 


( 4 ) 

( 5 ) 

( 6 ) 


where p is the density, v the velocity, P the pressure, $ 
the gravitational potential of the gas, and $sp the gravi¬ 
tational potential associated with the sink particles. For 
simplicity, we adopt an isothermal equation of state for 
the gas: 


P = CsP, 


( 7 ) 


where Cs = \/kT /(2.3m//) = 0.19km/s (P/lOK)^/^ is 
the isothermal sound speed of molecular gas with solar 
metallicity; T the temperature; and mn the mass of the 

hydrogen atom. _ 

The sink particle implementation of |Gong &: Ostriker] 
(2013) is briefly summarized here. Our sink particle cre¬ 
ation criteria require a local potential minimum and the 
density threshold pthr set by the Larson-Penston profile 
in Equation ([^ at a distance Aa:/2 from the center: 


Pthr = PLp(0.5Aa:) = 


8.86 


( 8 ) 


TT GAx' 

here Ax is the simulati on cell size. Previously, §4.3 of 


sink particle creation 

[20M 


et al. 


Federrath 


|Bate et al. 

19951 

al. 2(J1U Fac 

oan & 


Krumholz 

iNordlund 


20111), including the Truelove density threshold (Truelove 


et al.| 1997), additional checks for converging flow and a 


gravitationally bound state, and found no differences in 
sink particle evolution. In our algorithm, there is a mov¬ 
ing control volume (3Aa:)^ centered on the cell containing 
each sink particle that acts similar to the ghost zones of 
the simulation domain. At every time step, the density 
and momentum of the cells inside each control volume are 
reset by extrapolation of the surrounding non-sink zones. 
The gas mass and momentum accretion onto each sink 
particle are calculated based on the sum over all fluxes 
returned by the Riemann solver at the outer surfaces of 
the sink control volume. Flows onto sink particles due 
to shifting of the control volume are also included. Sink 
particle positions and velocities are advanced in time via 
a leapfrog integrator. The gravitational potential from 
and the forces on sink particles are comp uted using a 
particle mesh m ethod with a TSG kernal (| Hockney & 
Eastwood 1981). Sink particles merge if their control 
volumes overlap. 

2.1.2. Problem Specification and Units 

We used the prescription of convergin g super soni c flow 
with turbulent pe rturbations following |G011| and |Gong| 
& Ostriker (2013). The models represent a localized (~ 
pc scale) region m a turbulent GMG where the large scale 
velocity field converges and pre-stellar cores form in the 
post-shock dense gas. 

The model setup consists of supersonic flow augmented 
with turbulent velocity perturbations, converging to the 
mid-plane from +z and —z directions at an average veloc¬ 
ity of —Aics and PAics- The inflowing gas has the same 
spectrum of turbulent perturbations as the initial condi¬ 
tions (see below), i.e., the cloud “outer scales” are treated 
as maintaining a fixed turbulent amplitude. There is no 
explicit turbulent driving within the box, although tur¬ 
bulence can be driven by shock instabilities. In our pa¬ 
rameter survey, the Mach number Al values are 2, 4, 8 
and 16. The initial density is set equal to the code unit po 
(see below), as is the density of the converging flow from 
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the z boundaries. The boundaries in x and y directions 
are periodic. 

We apply velocity perturbations for both the whole 
domain initially and the i nflowin g gas subsequently, fol¬ 
lowing the prescription of |G01l| The scaling law of the 
velocity field represents the observed linewidth-size rela¬ 
tion i n GMCs over a range of scaled (|McKee & Ostriker 
20^ : 

5v{l) oc (9) 


However, the timescale for the turbulence to relax at the 
largest scales, attaining self-consistent density and veloc¬ 
ity perturbations, is longer than the inflow time from the 
boundaries. Thus, although the turbulence is not fully 
self-consistent, it serves to produce density perturbations 
that then evolve to collapse gravitationally. 

Assuming the cloud-scale supersonic turbulence pro¬ 
vides both the kinetic support of the cloud against its 
self-gravity and the large-s cale con verging flow with the 
virial parameter avir 1, GOll show (their Equation 
(43) and (44)) that the Mach number of the converg¬ 
ing flow can be related to the cloud’s size Rdoud and 
Jeans length Lj as M ^ Rdoud/Lj- With the scal¬ 
ing of the velocity in Equation the amplitude of ve¬ 
locity perturbations at the box scale L can be written 
as uid(T) = B This correspo nds to the 

“high-amplitude” perturbation simulations in |G01l| In 
this paper, we adopt 10% of the above values for veloc- 
ity per turbati ons, corresponding to the “low-amplitude” 
case in |G011[ Thus, the largest velocity perturbation in 
a simulation is limited by 


uid(T) =0.1(ML/Lj)1/2c,, (10) 

where L is the horizontal (x and y direction) box size. 
Note that for the Mach numbers and box size we use, al¬ 
though the average inflow velocity is highly supersonic, 
the velocity perturbations within the domain are sub¬ 
sonic. 

We choose to focus on the low-amplitude perturbation 
case for two reasons: Firstly, we wish to concentrate on 
the local details of collapse produced by converging flows. 
For realistic GMC conditions, A4 ^ 1 and the velocity 
perturbations at scales ~ Lj are small compared to the 
inflow motion Aics, which means a post-shock layer with 
simple geometry would form; this situation is amenable 
to local simulations. However, for completeness, we wish 
to explore a range of M. values down to A1 = 2, and at 
low M. the velocity perturbations at scales ^ Lj would 
be comparable to inflow speeds. In this case, there would 
no longer be a post-shock layer with simple geometry, so 
that global simulations would be more suitable than lo¬ 
cal. By reducing the perturbation amplitude, we can 
use local models for all Al. Secondly, we are limited 
by computational power in practice, and the GPU time 
required for a simulation run is increased in high am¬ 
plitude cases. Also, bigger simulation boxes and higher 
resolutions are required when the post-shock layer does 
not have a simple geometry and shows more complicated 
dynamics, increasing the necessary number of zones. 


^ Note, however, that at small scales where the turbulence is 
subsonic, it may follow a shallower Kolmogorov power law. 

^ In a global model, for which L —> 2itcioud, one can show that 
this reduces to EiD(2ftcioud) ~ Mcs- 


Nonetheless, it is important to address how the level of 
turbulence affects the core formation process. We there¬ 
fore conducted two sim ulatio ns with high-amplitude per¬ 
turbations (see Section 2.2 for model parameters), and 
compare to the low-ampiitude cases. 

The code unit of density po represents the average am¬ 
bient density in a GMC flowing into the region in which 
stars will form. For this reference density, a characteris¬ 
tic spatial scale for gravitational instability is the Jeans 
length 

1/2 


Li = c, 


( — ) 
\GpoJ 


The corresponding Jeans mass is 
Mj = poL^j = c] 




\G^PoJ 


and Jeans time is 




( 11 ) 


( 12 ) 


(13) 


We use the Jeans length, mass, and time at the unper¬ 
turbed density as the code units for length, mass, and 
time: Lq = Lj, Mq = Mj and to = tj- In making com¬ 
parison to observations, another useful quantity is the 
surface density integrated along the z-direction through 
the domain 

E = y p(a:,?/,z)dz = Soy (14) 

where Eg = PqLj. 

To convert code units of length, mass and time to 
physical units, one may choose an appropriate mean 
density, po = l-dn^f^ow// (where uh = ‘ 2 nH 2 ), nnd 
a temperature T. This allows consideration of GMCs 
with a range of properties. Many GMCs are observed 
to be self-gravitating, such that the virial parameter 
avir = ^RoMCcrldD/iGMcMc) is order-unity, where 
av,iD is the large-scale velocity dispersion along any 
direction. Defining the density of the cloud as po = 
Mgmc/{^'^Rgmc) 11^® surface density as Eqmc = 
4po.Rgmc/3 so that Rgmc = 5cr^pD/(avir7rGEGMc) , 
we can substitute pg = Jiravir G'^GMg/ (^Ocr^,!!)) and 
= -Mcs (assuming that the inflow Mach number 
M. is related to the large-scale GMC’s turbulence level), 
to obtain 


Ln — 


20 

3avii 


1/2 


M 


GE 


GMC 


Mg = TT 


/ 20 


y 3av 


1/2 


M 


G^Egmc’ 


(15) 

(16) 


and 


to 


20 

Sttvir 



GEgmc ’ 


(17) 


/ 3avir 
^ 20 


1/2 


1 

M 


Sgmc- 


(18) 


Eg = TT 
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These forms are useful because observations of GMCs 
measure surface density and Mach number more di¬ 
rectly than volume density. Many of the best-studied 
local dark clouds and inner-galaxy GMCs have surface 
density Sgmc ~ 2OOM0/pc , although clouds in more 
extreme environments such as the Galactic center and 
ultra-lu minous infrared ga laxies (ULIRGs) have higher 
Sgmc (|Dobbs et al.||2014|). Table [T] lists the two differ¬ 
ent ways of translating code units to physical parameters 
in molecular clouds. 

2.1.3. Core Finding Method 

The cores in the simulations are ident ified u sing the 
GRID core finding method developed by |GOlfj the re¬ 
gion belonging to a core is defined by the largest closed 
contour of the gravitational potential held around the 
corresponding local potential minimum that contains no 
other local minimum. Within each core, we further dehne 
the bound core region as all the material that has the sum 
of gravitational energy and thermal energy negative]^ 
The cores are often not spherical, so we dehne an effective 
radius of a core with volume Wore as Tcore = (^Wore)^^^- 
For the core properties such as mass, radius, and density, 
we use the notation Moore> ''’core) and Pcore for cores de- 
hned by gravitational potential alone, and Mcoreb, rcoreb) 
and Pcoreb for bound cores with Eth + Eq < 0. 

To compare the physical properties of cores at different 
stages of evolution, two sets of cores are identihed using 
different strategies: (1) In each simulation for a snap¬ 
shot at time ti when the hrst core collapses, the cores 
are identihed around each lo cal min imum of the gravi¬ 
tational potential, similar to GOll (2) For each sink 
particle that forms in the simulation, a corresponding 
core is identihed at tcoU, the time just before sink par¬ 
ticle formation. The cores in (2) are by dehnition grav¬ 
itationally unstable, and represent the collapsed stage 
when the asymptotic LP density prohle has developed, 
just before a protostar would form. These cores identi¬ 
hed by method (1) represent the structures in the dense 
regions of GMGs at an earlier stage of their individual 
evolution than the cores in method (2). Not all of the 
method (1) cores are bound, nor are they all guaranteed 
to eventually collapse. For brevity, we use the terms ti- 
core and tcoii-core respectively hereafter, for cores found 
via the diherent methods. 


2.2. Model Parameters 

To investigate any possible dependence on purely nu¬ 
merical parameters, we carried out a series of simula¬ 
tions with different box size and resolution for each Mach 
number. For each model with given box size and resolu¬ 
tion, we ran between 2 and 8 simulations with different 
random initial velocity held for the turbulence perturba¬ 
tions, so that statistics could be obtained for more than 
~ 100 cores. Because of the planer converging how ge¬ 
ometry, we chose our simulation box with = Ly > 

^ The gravitational and thermal energy density in each zone is 
Eq = —p{<I>max — 4*), where i&max is the gravitational potential 
at the core boundary, and Fth = 3l2pcl. We have also defined 
kinetically bound cores with the volume integral of Etj, + Ek + 
Eq less than zero, where Efc = [v^ F + v'^) 12 is the kinetic 
energy of the gas. We found that most cores are mainly supported 
by thermal pressure with > 1, even in tests with high- 

amplitude turbulence perturbation. 


to focus on the thin post-shock layer. For all simulations, 
Lz = 0.625Lj, and we have tested that using larger 
does not change the results significantly. The number of 
zones in each direction is set such that each cell is cubic: 
Lx/^x = Ly/Ny = Lz/Nz- The simulations were run un¬ 
til time tiim when the core formation rate starts to drop 
signihcantly because the material in the post-shock layer 
is most ly used u p by sink particle accretion. As shown in 
|G011| and also §3.1| in this paper, the cores form faster in 
higher Mach number simulations, and therefore we chose 
shorter tn^ for higher Mach number. 

Table summarises the model parameters for our 
study. 'Vtc mark in boldface the simulation set at each 
Mach number with the highest resolution and big gest 
box size, which will be used for further analysis in §3.1| 
and §3.3| if not specified otherwise. Models with smaller 
box size and resolution are mainly used to conhrm con¬ 
vergence, as discussed in p.2[ 


As discussed in Section 
have turbulence levels t 


roi 

at di 


the models of Table 
rat depend on Mach num- 
er following Equation (llOl). In addition, we have con¬ 
ducted two simulations analogous to model M08L2N256, 
but with turbulent perturbations an order of magnitude 
higher. 


3. RESULTS 
3.1. Overall evolution 
3.1.1. Development of Cores and Filaments 

First, we provide an overview of the dynamical evo¬ 
lution in our models. Figure shows the surface den¬ 
sity evolution of different Mach numbers at (2/3)ti, 
ti, (4/3)ti, fiini) where ti is the time for the first 
core to collapse in each simulation, corresponding to 

0.28to,0.23to,0.21to,0.13to for M = 2,4,8,16. The de¬ 
pendence of ti on M is discussed below. 

The notable density structures seen in all models are 
filaments and cores. Also shown in Figure [^are the loca¬ 
tions of sink particles. Because gaseous rilaments are 
unstable to long itudinal self-gravitating fragmentation 
(Nagasawa 1987), the suggestion has been that there is 
a two-stage process for core formation, with hlaments 
first forming, and t hen fragmenting gravitationally (e.g., 
Andre et al. 2014). However, we find a subtly differ¬ 
ent progression! Figure makes clear that cores and 
filaments develop simultaneously, rather than hlaments 
forming hrst, followed by fragmentation into cores. Snap¬ 
shots at closely-space time intervals show that the central 
densities of proto-cores and hlaments grow together, un¬ 
til some of the gravitationally unstable cores collapse and 
form sink particles. As a result, the sink particles are not 
randomly distributed in the post-shock region, but are 
dotted along the long thin hlamentary over-dense struc¬ 
tures. Because our simulations do not include feedback 
to limit accretion, sink particles continually accrete sur¬ 
rounding material, and after some time may merge with 
one another. 

The similarity in the density structure patterns of Ml = 
2,4 or Ad = 8,16 arises because they have the same 
seeds for turbulence perturbations. Because of the planar 
converging how geometry, the x and y component of the 
large scale velocity are not changed by the shock, and 
these velocities seed the density huctuation patterns that 
grow to create hlaments and cores. Thus, the structure 
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TABLE 1 

Code units and corresponding physical units 
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Note. — Code units are given as a function of reference volume density njf Q (left) and based on a corresponding cloud surface density 
Eqmc Mach number Ai (right). 


TABLE 2 

Simulation Set Model Parameters 


^Name 

M 

number of simulations 

box size (I/q) 

resolution 

flim (^o) 

M02L6N256 

2 

3 

6 X 6 X 0.625 

1536 X 1536 X 160 

0.5 

M02L4N256 

2 

6 

4 X 4 X 0.625 

1024 X 1024 X 160 

0.5 

M02L4N128 

2 

4 

4 X 4 X 0.625 

512 X 512 X 80 

0.5 

M02L4N64 

2 

8 

4 X 4 X 0.625 

256 X 256 X 40 

0.5 

M04L6JN256 

4 

2 

6 X 6 X 0.625 

1536 X 1536 X 160 

0.4 

M04L4N256 

4 

3 

4 X 4 X 0.625 

1024 X 1024 X 160 

0.4 

M04L4N128 

4 

3 

4 X 4 X 0.625 

512 X 512 X 80 

0.4 

M08L2N512 

8 

8 

2 X 2 X 0.625 

1024 X 1024 X 320 

0.35 

M08L2N256 

8 

6 

2 X 2 X 0.625 

512 X 512 X 160 

0.35 

M08L2N128 

8 

5 

2 X 2 X 0.625 

256 X 256 X 80 

0.35 

M08L4N256 

8 

2 

4 X 4 X 0.625 

1024 X 1024 X 320 

0.35 

M16L2N512 

16 

4 

2 X 2 X 0.625 

1024 X 1024 X 320 

0.2 

M16L1.5N512 

16 

6 

1.5 X 1.5 X 0.625 

768 X 768 X 320 

0.2 

M16L2N256 

16 

3 

2 X 2 X 0.625 

512 X 512 X 160 

0.2 


^The naming convention is based on the Mach number, the horizontal box size in units of Lg, and the number of zones in Lq. 


provided by turbulence is crucial to initiation of the self- 
gravitating cores that collapse to make sinks. However, 
it is not true that the cores are formed simply from the 
growth of initial perturbations. Filaments and clumps 
can merge and fragment before they finally collapse into 
sink particles. 

For the present simulations, the initial turbulent mo¬ 
tions are generally subsonic. Although the structures 
that grow are seeded by the turbulent compression, self¬ 
gravity is crucial to their developmenliN This is true for 
both filaments and cores. We note, however, that on 
larger scales for which flow velocity would be supersonic, 
filamentary structures can form from turbulent compres¬ 
sion alone. This sort of filament, for which gravity is not 
important, is distinct from the filaments formed in our 
simulation^ 

Figure shows an x = const slice through the density 
maximum at the time when collapse leads to the first sink 
particle. Evidently, the post-shock layer seen in cross 
section can be highly nonuniform, especially for Ad = 8 
and 16. This is due t o nonlinear thin -shell instabilities of 
the post-shock layer ( Vishniac|1994 ) that are stronger at 
higher Mach number. As a result of these instabilities, 
the gas is vertically dispersed and the density in the post¬ 
shock region is lower than the value that would 

apply for a simple isothermal shock. The thickness of 
the post-shock layer is correspondingly higher than the 
value that would apply in the absence of instabilities. In 
fact, at later times, the median density (and pressure) 
of the post-shock region is nearly independent of Mach 


^ Even for the models with larger-amplitude turbulence, gravity 
is i mportant for creatin g core-forming filaments. 

® |Andre et al.| ||2010| discusses observations of filaments in both 
non-selt-gravitatmg and self-gravitating clouds. 


number. This leads to a very shallow dependence of core 
mass on Mach number, which is discussed further in §3.3[ 
We note that mag netic fields reduce insta bilities of the 
post-shock region (Chen & Ostriker 2014), so that the 
post-shock density and magnetic pressure increase with 
inflow Mach number. 

The map of Vx at different times is shown in Figure 
At the early core-building stage, the velocities are 
sub-sonic, con verging towards the ridge of g rowing hla- 
ments (see also G01l[ Chen & Ostriker 20141, due to the 


gravity from cores and hlaments. Alter a core collapses 
to reach an LP profile, a sink particle is created and the 
surrounding material accretes onto the sink particle at 
the free-fall velocity. 


3.1.2. Collapse of the First Cores: the Non-linear Time and 

Mass 

Treating the post-shock layer as a slab with half thick¬ 
ness H and surface density E = 2csMpot that grows 
in time, f or linea r-amplitude in-plane perturbations with 
kH 1, GOll showed the wavelength of the unstable 
mode that would have the greatest exponential amplifi¬ 
cation at a given time is 


A™ _ / 2a/3 \ 1 

Lo Ad 1/2’ 


(19) 


where F^ax = lii(5Eniax/(5Einit- The time at which an 
amplification Fmax is reached is 


tm 

to 


/2F„ 


V 


1/2 


(20) 


Adi/2' 
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Fig. 1. — E volution of surface density projected along the direction of converging flow (£). Logarithmic color scale shows S/Eq (see 
Equation i Mi and Table li for Mach number 2, 4, 8, 16 at times 2/3£i, ti, 4/3ti, fiim- The x- and y- axes are in units of Lq = Lj (see 
Equation iTTi and Table Ti. For A4 = 2,4, we show only a part of the simulation domain (a square region with a size of 2Lo) for direct 
comparison with A4 = 8,^ simulations. A4 = 2,4 and A4 = 8, 16 simulations have the same seeds for turbulent perturbations, and this 
leads to similar large-scale structure. The sink particles are over-plotted as spheres with color scales indicating their masses (note that the 
plotted size exceeds the size of the sink control volume). 


A characteristic mass associated with this mode is 
_ ( 2^3 1 

Mo Mo ^ ^ 

We fitted the time for the first core to collapse ti 
and the average mass of the first five cores in each 
simulation (Mgi-stcores) as a function of Mach number. 
In addition to the primary models (boldface in Table 
we also included M02L4N1024, M04L4N1024, and 
M16L1.5N768 models, which have the same resolution 
and have been demonstrated to have sufficient box-size 
(see |3.2| ). The result is shown in Figure]^ Our fit gives 
ti/to = and (Mfi 

rstcores )/Mo = 0.38M-°'63. 

The dependence of ti and (Mg^stcores) on Mach number is 
similar to the dependence tm, M^ oc of Equations 


(20) and (21). and Fmax = 1 — 2 would give coefficients 
in Equations (20) and (21) comparable to our fits. 

However, theiarge dispersion in (Mfij-stcores) (see Fig¬ 
ure]^ suggests a process more complex than local growth 
of cores with in-plane waves lx ^ ly ^ Am- In particular, 
filament development with lx and ly unequal clearly plays 
a role in core formation. We note that using Fmax = 1) 
gives Am/To = 1-86/Al^^^, which we find is 


19 


Equation 

roughly equal to the typical filament separation se en in 
Figure [T| and also consistent with what was found in|Van| 


Loo et al. (2014). Thus, both the time to first collapse 
and the prominence of filaments implies that non-linear 
instability of asymmetric self-gravitating in-plane modes 
is important to core formation. 

Finally, we note that the mass per unit length associ¬ 
ated with the most amplified mode is AmS(tm) = dCg/G, 








































Fig. 2. — Slice of density along the x-axis at the time and position of the first sink particle. The sink particle in each panel is marked with 
a black dot. Notice that the higher A4 models have strongly dispersed post-shock layers (along the inflow direction i), as a consequence of 
hydrodynamic instability of the shock-bounded layer. 


which is equal to twice th e critical mass per unit length 
for an isothermal filament ( Ostriker]l964 ), or 33Mqpc~^ 
for T — IOK]^ Of course, this mass is not all available at 
time tm, because only a fraction of this has been gathered 
into the filament. We find that typical values of the mass 
per unit length in filaments at time ti is ~ (1 — 3)Cs/G. 

3.1.3. Evolution of Cores 

Examples of core finding for individual tcoii-cores and 
an overall inap of ^1 ti-cores in one simulation are shown 
_ We can see clearly that cores are 
aments. A typical mass core is often 
found embedded in dense hlaments. Although high mass 
cores may form in relative isolation early on, these loca¬ 
tions often become junctions of filaments later on. These 
structures of cores lying in filaments or at the junctions 
are qualitatively very similar to the observation of star¬ 
forming regions with He rschel (lAndre et al.|2014 |, as well 
as previous work (e.g., Johnstone & Bally 1999 Hart- 
maim||2002 ). 

|GUU9| classified the core development into four differ¬ 
ent stages: core building, core collapse, envelope infall, 
and late accretion. These stages can also be seen in our 
simulations, as illustrated in Figure[^and[^of the density 
(angle-averaged) and velocity promes of^typical cores. 
Initially at the core building stage, the core evolves in 
quasi hydrostatic equilibrium, with subsonic internal ve¬ 
locities. The core mass and density grow slowly as the 
gas from surrounding environment flows into the core 
potential well. When the core becomes gravitationally 
unstable, it collapses in a short timescale, reaching the 
LP profile. Then a sink particle is created in the core 
center, and starts accreting the gas from the envelope, 
causing the density to drop as the core collapses inside- 
out and approaches p(r) oc r~^-^ a proHle close to the 
expectation for free-fall ( Shu||1977 ). 

Noticeably, the angle-aver aged core density profiles 
continue to drop smoothly beyond the effective core ra¬ 


dius Tcore = [3I4ore/(47r)] as shown in Figure and 
We note, however, this is partly due to the method 
we use to calculate p(r), which takes the average density 
in a spherical shell at distance r from the core center, 
including the low-density pre-shock gas at t/Lq > 0.04 
outside of the planer post-shock layer. For example, the 
density profile of the core in Figure only extends to 
^ 2rcore along some directions and is much less smooth 
than the angle-averaged density prohle, as can be seen 
in the upper right of Figure]^ The anisotropy of cores 
(e.g., as in the upper left of Figure can also make 
the angle-averaged density differ from the density profile 
along individual principal axes. Nevertheless, cores often 
extend beyond their effective radius Tcore; because they 
are formed in denser environments than the average post¬ 
shock density, and instantaneous tidal forces (which de¬ 
fine Tcore) do not limit the material that may ultimately 
fall into a sink particle. In fact, we often see that the 
mass of the sink continues to grow after it reaches the 
mass of the tcoU-core. There is no obvious boundary be¬ 
tween the en velope infall and late accretion stages, as 
discussed in §3.4.1 

In our tests of the case with high-amplitude pertur¬ 
bations, the overall evolution is very similar to the low- 
amplitude cases. Cores and filaments still form in a sim¬ 
ilar manner, with cores embedded within filaments. At 
smaller scale {I LA, both the initial velocity pertur¬ 
bation in Equation is smaller and the turbulence dis¬ 
sipation timescale t = l/v (x is shorter. As a result, 
even when the large-scale velocity perturbation is super¬ 
sonic, the high density regions within hlaments where 
cores form are generally sub-sonic or trans-sonic, lead¬ 
ing to very similar core evolution dominated by thermal 
pressure and gravity. Individual cores formed in simu¬ 
lations with high-amplitude velocity perturbations show 
the same internal collapse and infall prohles as shown in 
Figures and 

3.2. Convergence of isothermal fragmentation 


® For a self-gravitating isothermal sheet with fixed surface den¬ 
sity S and scale height ("/(ttGS), the mass per unit length of the 
fastest growing mode is also ^ Jc^/G. 


In order to understand protostar formation in our 
simulations, it is important to test numerical conver¬ 
gence, both to make sure the fragmentation process is 
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Fig. 3.— Color scale shows maps of Vx averaged over the z-axis; (vx) = f Vxpdzj f pdz, for one j\4 = 8 simulation at t/ti = 
2/3,11/12,1,4/3. The white and magenta curves are contours of surface density S/Sq = 8 and 24. Sink particles are marked as black 
dots. Gas converges toward the ridges of forming filaments, and then accretes at free-fall onto sink particles. 




Fig. 4.— Time of first collapse ti (left) and average mass Mfirstcores (right) of the first cores verses Mach number j\4. The dots are 
the average values for all the included models (see text), with error bars showing the standard deviations. The dashed line is the log 
linear fit, and the dotted line is t he f it wit h a fixed slope of -1/2. The fitting gives ti oc and (Mfirstcores) cx; sim ilar 

to tm,Mm oc in equation l |20[ | and ( [2l]l . If using a fixed slope of -1/2, the fitting gives Fmax = 1.3 for ti in equation ( |20[ | and 

Fmax = 2.4 for (Mfirstcores) equation ( |21[ l. ine green solid line in the right panel plots Mm in equation ( |21[ l with Fmax = 1* 
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M=2 


< 0.0 




x/Lq 


Fig. 5.— Illustration of individual tcoircores, for A4 = 2 (left 
panels) and A4 = 16 (right panels) simulations. The upper and 
lower panels show typical-mass and high-mass cores res pect ively, 
which are marked as red dots and star symbols in Figure |12| The 
color scheme shows surface density log;i^Q(S/So). The black and 
magenta curves draw the boundaries projected along the z-axis 
of cores defined by gravitational potential alone, and bound core 
defined by Eth + 0 (not plotted in the upper panels for 

clarity). 



- 0.55 


< 0.30 


Fig. 6.— Illustration of ti-cores for one A4 = 8 simulation. Black 
and magenta curves are as in Figure]^ 


not purely numerical, and the cores in our simulations 
are well resolved. In the literature, there has been a 
debate of whether fragmentation in isothermal HD sim¬ 
ulations is determined by pnrely numerical rather than 
physical factors. If an isothermal sphere were to col¬ 
lapse homologously, the Jeans length within the sphere 


would decrease with radius as oc R' 


3/2 
sphere' 


This led to 


the early idea ( Hoyle|1953 1 of hierarchical fragmentation 
into smaller pieces until pieces are no longer isothermal. 
In reality, however, core collapse is highly non uniform 


and leads to centrahy-concentrated structures (Boden- 


heimer & Sweigart 

1968 

Larson 

1969 

Pension 

1969 


length is comparable to the radius. Thus it is not obvi¬ 
ous that fragmentation can happen at arbitrarily small 
scales even for a simple isothermal equation of state. 


Martel et al. (20061 carried out a series of isothermal 
SFH simulations employing particle splitting techniques, 
and argued that the core masses they found are deter¬ 
mined by numerical resolution and the density thresh¬ 
old of sink particles. However, the highest resolution 
they used only marginally resolves the Jeans mass at the 
density threshold for their sink particles, and they did 
not test with a higher resolution or density threshold. 


Krumholz (20141 further suggested that isothermal evo- 
lution with selt-gravity cannot lead to a characteristic 
distribution of fragment mases, and that stellar masses 
must depend crucially on non-iso thermal effects. How¬ 
ever, Inutsuka & Miyama (1997) found that filaments 
(even with supercritical mass per unit length > 2c^/C?) 
fragment longitudinally (rather than collapsing to a spin¬ 
dle) provided that the initial density perturbations are 
sufhciently large. In real clouds, as cores and filaments 
begin their growth together, filaments may not become 
supercritical until after cores have become nonlinear. 

We find that at the stage of core collapse, there is in 
fact a well defined characteristic mass scale for the ideal 
isothermal case. Figure shows the core mass function 
(CMF) for tcou-cores of models with different resolutions 
for each Mach number. For the lower resolution models, 
the distribution of core mass drops off rapidly at rcore ^ 
2Ax. This shows that cores with < 12 zones may not 
be well resolved. However, as we increase the resolution, 
although the CMF continues to extend towards lower 
masses, thepeak of the CMF does not vary much (see 
also Table [3 ). We note that the density threshold for 
sink particle creation also increases with resolution as 
described in Equation (1^. For the highest resolution 
model of each Mach nurnber, the CMF peaks at a mass 
of cores that are well resolved with rcore/Ax = 7—13. 
Although simulations with even higher resolutions could 
show a modification of the CMF, especially at the low 
mass end, our models suggest that the CMF peak, or the 
characteristic mass of the core, would remain the same. 

We also investigated the numerical effect of horizontal 
box size L in our simulation. The box size limits the 
longest wavelength perturbation that can grow in the 
simulation domain X < L, and also limits the amplitude 
of initial vel ocit y perturbations as described in Equation 

S . Figure[lo|shows histograms of Mcore of models with 
;rent box sizes. The distributions of Mcore do not vary 
appreciably with increasing box size in our models (see 
also Table . This shows that for our simulations with 
a box sufnaently larger than the do minant non-linear 
wavelength (L > Am, see Equation 19), the box size does 
not have a significant effect on core formation. 

Although our simulations show a well-defined, con¬ 
verged value for the peak in the distribution of core 
masses at the point of singularity formation, fragmen¬ 
tation at later stages could in principle further alter this 
distribution. In particular, infall of rotating envelopes 
leads to formation of prestellar accretion disks, which 
we do not model in the current study. These accretion 
disks would be susceptible to fragmentation if external 
and internal heating are not properly taken into account. 
We conclude that it is important to distinguish between 
the core stage and subsequent stages in assessing the 
outcome of self-gravitating fragmentation in turbulent 
clouds. Going from the CMF to IMF physically involves 
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Fig. 7. — D ensity and velocity evolution of a typical mass core in one A4 = 2 simulation (marked as a red sphere in the upper left panel 
of Figure [T^ and shown in the upper left of Figure]^. The blue, cyan, green, orange and red lines show the core profiles at t/tg =0.260, 
0.343, 0.3Y3 (tcoii)i 0.403, and 0.467. The crosses mark on each profile the grid centers in the simulations. This sink particle will merge 
with another sink at fmerge/to = 0.483. Upper left: angle-averaged density along r. The black dashed line shows the LP profile (Equation 
[^. The red dashed line is a log-linear fit of the red density profile for 2Ax < r < 2rcore, which gives p(r) oc y—1-56 Upper right, lower left 
and right: v^, Vy and Vz along the x, y and z axis. The orange and red dash-dotted lines show the free-fall velocity Vff = y^2GMap/r, 
where Map is the mass of the sink particle at the corresponding time of the orange and red velocity profiles. The dotted vertical lines 
denote the value of rcore = [3V'core/(4'7r)]^/® at tcoll- 


structures at increasingly high temperature, and simula¬ 
tions must include these non-isothermal effects to model 
the stages of disk formation and evolution. 

3.3. Core Properties 

3.3.1. Typical Core Mass, Radius, and Density 

We summarise the basic physical properties for tcoir 
cores and ti-cores of different models in Table and 
Table Cores containing fewer than 27 grid cells are 
considered to be poorly resolved (similar to the resolu¬ 
tion limit in Figure!^ and are not included in Table 
or further analysis. Tableagain shows the convergence 
of the characteristic mass for simulations with different 
resolutions and box sizes, as discussed in §3.2[ 

The median tcoii-core sizes and masses are quite simi¬ 
lar for different Mach numbers, as a result of the sim¬ 
ilar post-shock density across different Mach numbers 
Ppost-shock/po 10 — 20. Using the fiducial values of 
Mq and Lq from Table [B the characteristic core masses 
and radii are ~ (1 — 2)Mq and ~ (0.02 — 0.03)pc. Fig¬ 
ure [m plots the mass defined by gravitational potential 
alone of tcoii-cores (Mcore) versus the post shock criti¬ 
cal Bonner-Ebert sphere mass (Mbe) a.nd Mach number 


(M). Setting Pedge = c^/Opost-shock in Equation ([^, the 
critical Bonner-Ebert mass becomes 

Mbe — q 22 f —shock 

Mq \ Po 



We fit the post-shock density field (defined as zones with 
p > l.Spo) with a log-normal distribution with mean 
value of log]^g(p/po) equal to p and dispersion a. This 
distribution remains roughly constant over the time of 
core formation. We then used p, and a at the time when 
the core formation rate is roughly at its peak to estimate 
the corresponding distribution of Mbe'- from Equation 
lUBB(ppost-shock)/Mo also follows a log-normal dis¬ 
tribution with a mean value log]^Q(1.2/7r^/^) — 0.5/i and 
dispersion O.Str, which is plotted as dots and error bars 
on the x-axis in Figure 11 On the y-axis of Figure [TH we 
plot the median value of Mcore a.nd error bar showing the 
median absolute deviation of Mcore in logarithm space, 
as listed in Table Typically 0.5 < Mcore/M be 2, 

consistent with the general range expected for isothermal 
fragmentation. The median core mass depends on inflow 
Mach number as Mcore ~ As noted above, the 
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Fig. 8. — Same as Figure|^for a typical core in one A4 = 16 simulation run (marked as a red sphere in the lower right panel of Figure [T^ 
and shown in the upper rigm of Figure]^ at t/to =0.100, 0.134, 0.146 (tcoll)i 0.151 and 0.165. This sink particle will merge with another 
sink at imerge/to = 0.171. Upper left: The red dashed line gives a fit of p(r) oc y—for the density profile during the infall stage. The 
velocities are more turbulent than in the j\4 = 2 model, but still mostly subsonic within the core before it collapses. 

weak variation in Ppost-shock, and resulting weak depen¬ 
dence of the characteristic mass on Af, is a result of 
instability of the post-shock layer. In magnetized simu¬ 
lations, for which these instabilities are suppressed, the 
post-shock density is close to the expectation of magne¬ 
tized isothermal shock Ppost-shock tx A4. Combining this 
effect with anisotropic flows along magnetic field lines 
leads to a stronger dependence of core mass on Mach 
number close to Mcore oc (Chen & Ostriker 2015, 

in preparation). 

With increasing Mach number, we find that more cores 
form per unit area in each simulation run. Our fitting 
gives ncoie/{L/Lo)‘^ ^ for tcoircores, similar to 

what would be expected from the decreasing separation 
of filaments oc of Equation (19). 


In comparison to the tcoircores, the -cores show a 
wider range of physical properties. This can be explained 
by the fact that they include cores at very different evo¬ 
lutionary stages, as discussed in detail in the following 
sections. The median masses of cores are lower, and me¬ 
dian radii are higher, for ti-cores compared to tcoii cores. 
The properties of ti-cores also show a stronger depen¬ 
dence on Mach number, because they represent the non¬ 
linear size and mass in t he p ost-shock sheet that scales 
as (see Equat ion ( 19) and (21)), as also shown in 

the earlier studies bylCUllJ 


Figure 12 plots the mass and time of collapse for ev¬ 
ery tcoircore. For higher Mach numbers, the cores form 
so oner an d in a less extended period of time, as explained 
in §3.1. 2| This also leads to more similarity between ti 
and tcoii cores (see Table and for higher Mach num¬ 
bers. There is a similar pattern m models with different 
A4. The average core mass decreases slightly with time, 
as smaller cores form in dense filaments which become 
more developed at later times. At the end of each simu¬ 
lation, the core formation rate drops, when most material 
in the post-shock region has been accreted by the sink 
particles. This is further discussed in §3.4.2[ 

The median value and the median absolute devi¬ 
ation at a logarithm scale of core mass and radius 
in high-amplitude velocity perturbation simulations are 
logio(Mcore/Alo) = -1.29 ±0.22 and iog^^ir^oTe/Lo) = 
— 1.86 ± 0.18. This is about ~ 60% of the core mass 
and radius in low-amplitude perturbation cases (see the 
M08L2N256 model in Table|^, but subject to the limited 
statistics of 21 cores. There are also a few massive tur¬ 
bulence supported cores formed in high-amplitude simu¬ 
lations, which are not seen in low-amplitude cases; this 
is further discussed in the last paragraph of §3.3. 4| 

3.3.2. Core Binding 

Since the tcoii-cores are identified at the time of collapse 
by definition, they are expected to have density profiles 
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Fig. 9. — The core mass function (CMF) for fcoircores of models with different resolutions at each Mach number. Here the cores are 
defined by gravitational po tential alone. The y-axis is normalized such that the area of each histogram is one. The vertical lines are 
Mi^p{2Ax) in equation ( |23[ ), on the left of which the cores are considered not well resolved. There are about 140-300 cores in each model. 
Each histogram is normalized so that the total mass equals Mq. Although the minimum core mass decreases with increasing resolution, 
the peak of the CMF does not. 


TABLE 3 

Summary of Physical Properties for tcoii-coRES. 


Model 

^fo§10 (-^core/Afg) 

Lq) 

^fologio (pcore/Po) 

^core,total 

'^ricore/(L/Lo)^ 

M02L6N256 

-1.07 ± 0.29 

-1.47 ± 0.33 

2.70 ±0.66 

267 

2.47 

M02L4N256 

-1.08 ±0.27 

-1.45 ±0.28 

2.63 ±0.55 

239 

2.49 

M02L4N128 

-0.99 ±0.17 

-1.41 ± 0.21 

2.64 ± 0.44 

169 

2.64 

M02L4N64 

-0.82 ± 0.14 

-1.23 ±0.14 

2.27 ± 0.29 

220 

1.72 

M04L6JN256 

-1.09 ±0.22 

-1.42 ± 0.21 

2.56 ± 0.42 

346 

4.81 

M04L4N256 

-1.12 ± 0.21 

-1.46 ±0.21 

2.67 ± 0.45 

266 

5.54 

M04L4N128 

-1.08 ±0.15 

-1.41 ± 0.17 

2.58 ±0.36 

264 

5.50 

M08L2N512 

-1.18 ±0.28 

-1.55 ± 0.29 

2.84 ±0.60 

232 

7.25 

M08L2N256 

-1.17± 0.18 

-1.51 ± 0.19 

2.70 ± 0.42 

197 

8.21 

M08L2N128 

-1.15 ± 0.14 

-1.49 ±0.14 

2.70 ± 0.32 

182 

9.10 

M08L4N256 

-1.17 ±0.19 

-1.50 ±0.20 

2.70 ± 0.44 

121 

7.56 

M16L2N512 

-1.26 ±0.20 

-1.57± 0.21 

2.87 ± 0.48 

122 

7.62 

M16L1.5N512 

-1.24 ±0.23 

-1.53 ±0.25 

2.71 ± 0.51 

107 

7.93 

M16L2N256 

-1.33 ±0.17 

-1.62 ±0.19 

2.91 ± 0.42 

158 

13.17 


^Median value ± median absolute deviation (MAD) 

^ pcore = Mcore/Vcore for each core. 

^Average number of cores per simulation run per Lq area in the x — y plane. 
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Fig. 10.— Similar to Figure]^ for models with different box sizes. 


TABLE 4 

Summary of Physical Properties for ti-coRES. 


Model 

^®Sl0 i^core / Afo) 

loglo(’'coreAo) 

^Ogio(pcore/ po) 

^core.total 

^core/ 

M02L6JN256 

- 1.33 ±0.38 

-1.16 ± 0.15 

1.53 ± 0.09 

94 

0.87 

M02L4N256 

- 1.44 ± 0.35 

-1.24± 0.12 

1.59 ±0.14 

86 

0.90 

M04L6JN256 

-1.51 ± 0.24 

-1.31 ± 0.09 

1.83 ± 0.12 

151 

2.10 

M04L4N256 

-1.43 ±0.27 

-1.32 ±0.10 

1.82 ±0.16 

103 

2.15 

M08L2JN512 

-1.57 ± 0.32 

-1.44 ± 0.13 

2.00 ± 0.17 

139 

4.34 

M16L2N512 

-1.60 ±0.29 

-1.51 ± 0.18 

2.36 ± 0.29 

103 

6.44 

M16L1.5N512 

-1.46 ±0.28 

-1.48 ±0.18 

2.30 ±0.26 

83 

6.15 


Note. — Notations same as Table [s] 
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Fig. 11.— Left: Median mass of tcoircores (Mcore) compared to the post shock critical Bonner-Ebert sphere mass {Mbe) for different Mach 
numbers (see text for definitions). The dashed line shows Mcore = Mbe the shaded area marks the region with 0.5 < McorelMBE < 2. 
Right: Mcore versus Ai. The dashed line shows a fit of the median Mcore with Al, giving Mcore/Mo = 0.097A1“^’^^. 
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Fig. 12.— Mass and time of collapse of individual tcoircores for different Al. The red dots and stars mark a typical mass core and a 
high mass core in Ad = 2,16 simulations as illustrated in Figure]^ 


close to the LP profile. This is indeed the case, as shown 
below (see also Figure[^and|^. Integrating Equation ([^, 
the mass of the core ■^1 be proportional to its radius 

8.86 r 

Mq tt Lq ' 

Solving the Poisson equation 


(23) 


1 d 


,d$\ 




(24) 


the gravitational potential of an isolated core with LP 
profile can be written as 


- A$(r) = $edge - $(?’) = -8.86cf In ■ 


(25) 


where $edge is the gravitational potential at the edge of 
the core, r^oie- Using = — A<i)(r(.oreb) for the definition 
of bound cores, we can obtain the relationship of ^coreb 
and Tt-ore: 


^coreb — I’core^ 


- 1 / 8.86 


0.9rc 


(26) 
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The left panels of Figure shows the mass-radius re¬ 
lation of the tcoU-cores defined by gravitational poten¬ 
tial alone and the bound tcoii-cores with Eg + E th < 
0. The results are consistent with equation (231, i.e. 
Mcore OC Tcore and M^oveh oc Tcoreb, showing thaFessen- 
tially all cores approach the LP density profile at the 
time of collapse. The left panel of Figure plots the 
size of every tcoii-core defined by gravitational potential 
alone compared to the corresponding bound core with 
Eg + Eth < 0, showing that rcoreb oc ?'corej as expected 
from Equation (261. Due to the potential of surround¬ 
ing cores and filaments, gravitational potential profiles at 
the edges of the tcoii-cores in the simulations are slightly 
flatter than that they would be for isolated cores. There¬ 
fore, we Hnd rcoreb/^’core ~ 0.75, slightly smaller than the 
coefficient in Equation (26). 

In comparison, many ti -cores are not strongly cen¬ 
trally concentrated or gravitationally bound, as shown 
in the right panels of Figure and Figure The 
mass radius relation of ti-cores is found to be iW oc 
with fc = 1.2 — 2.5, implying that they have flatter den¬ 
sity profiles than the LP profile or critical Bonner-Ebert 
sphere. This is consistent with results for ti-cores i n pre¬ 


vious simulations with or without magne tic helds (Gong 
fc Ostrik^|2009 Chen fc Ostriker||2014 ), and similar to 


many core-property surveys i n diHerent molecu lar clouds 
(e.g., Curtis fc Richer 2010 Kirk et al. 2013), in which 
k = 1.4 — 2.4 with various tracers (see h'igur e 7 and 
corresponding discussions in Kirk et al. (2013|). This 
suggests that cores identified at any given time in sim¬ 
ulations or observations (C-cores in our case) can be in 
very different stages of evolution. Because any individ¬ 
ual core spends most of its lifetime in the core building 
stage, many of the cores in a given snapshot of time will 
resemble sub-critical Bonner-Ebert sphere^ with only a 
fraction of them more evolved and in the core collapse 
stage. Some of the cores we identify, especially low mass 
cores, are structures with weak self-gravity and may not 
finally collapse. The spread of mass-radius relation is 
larger for lower mass cores in the right panels of Figure 
[Tsl implying a wider range of evolutionary stages. Inter¬ 
estingly, we find a shallower mass-radius relation at ti for 
higher Mach number simulations, indicating those cores 
are more evolved and therefore more centrally concen¬ 
trated. With complete population studies from Herschel 
and ALMA that are able to quantify the relative popu¬ 
lations of strongly concentrated cores versus those with 
flatter density profiles, the timescales for core c ollapse 


versus core building stages can be measured (see Marsh 


et al. (2014) for initial Herschel results). 


Another way to quantify whether a core is gravitation¬ 
ally unstable or not is to look at the ratio of core mass to 
the critical Bonner-Ebert sphere with the same average 
density of the core. For the LP density profile, the ratio 
Mcore/MBsip = pcore) = 6.8 (see Equation ^). In Fig¬ 
ure all fcoii-cores are gravitationally unstable, with 
the ratio Mcore/MBsip = Pcore) > 1- In comparison, 
among all the ti-cores , only ~ 45% of cores defined by 
the gravitational potential alone, and ^ 60% of bound 
cores with Eg + Eth < 0, are gravitationally unstable 


^ Cores at very early building stages would, however, be difficult 
to pick out in a noisy background, as the density contrast would 
be low. 




Fig. 13.— Mass-radius relation of cores. Left panels: fcoir 
cores defined by gravitational potential alone (upper) and bound 
tcoircores with Eq + F^th < 0 (lower). The solid line plots th e 
mass-radius relation for the LP density profile in Equation ( |23[ l. 
Right panels: same as for left, but for ti-cores. Dotted lines plot 
fits of mass-radius relation M oc for different A4, and the values 
of fitted k are listed on the plot legends. The ti-cores are less 
internally stratified than tcoii-cores. 



Fig. 14.— Radius of cores defined by gravitational potential 
alone (rcore) compared to radius of the portion of the bound cores 
with Eq Eth < 0 (rcoreb)- Left: tcoii-cores. The dashed line 
gives a linear fit of rcoreb/^core = 0.75. Right: ti-cores. The solid 
line indicates Tcoreb/'^core = 0.75. Most cores at ti are not strongly 
bound. 

with Mcore/MBE{p = Pcore) > 1- There is also a clear 
trend that this ratio increases with mass, which is not 
present for the fcoii-cores. This again suggests that many 
ti-cores, especially the low mass ones {Mcore/AIo < 10“^ 
in the middle and right panels of Figure [T^, are transient 
structures from turbulence perturbations, and may not 
finally give rise to star formation. This is also evident in 
the CMF described in §3.3.4| 

The core binding properties in high-amplitude pertur¬ 
bation simulations are similar to the low-amplitude cases 
as discussed above. All of the tcoii-cores are strongly 
gravitationally bound, whereas most of the ti-cores, es¬ 
pecially the lower mass ones, are not strongly concen¬ 
trated or highly stratihed. 

3.3.3. Core Shape 
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Fig. 15.— The ratio of core mass to the mass of a critical Bonner-Ebert sphere with the same internal average density, versus core mass. 
Left:tt,oii-cores defined by the gravitational potential alone. Middle and right: ti-cores, defined by gravitational potential alone (middle), 
and bound cores with Eq + Eju < 0 (right). The magenta and black dashed lines are Mcore/MBEiP = pcoie) = 6.78 (for LP profile), and 
1, respectively. 



Fig. 16.— Distributions of three-dimensional core aspect ratios of 
cores defined by gravitatio nal pote ntial alone for each Mach num¬ 
ber. See also Figure 19 in |G011| Cores lying on c/a = b/a are 
formally prolate and along b/a = 1 are formally oblate. We sub¬ 
divide (see diagonal lines) and classify as follows: approximately 
prolate (between c/a = 1 and c/a = 1.5b/a —0.5), triaxial(between 
c/a = 1.5b/a — 0.5 and c/a = 3b/a — 2), and approximately 
oblate(between c/a = 3b/a — 2 and b/a = 1). Left: tcoll-cm’es. 
Approximately 66% of icoii-cores are prolate, 29% triaxial and 5% 
oblate, independent of Mach number. Right: ti-cores. Approx¬ 
imately 30 — 50% of ti-cores are prolate, 30 — 40% triaxial, and 
10 — 20% oblate. Slightly more cores are prolate for higher Mach 
numbers. 

The shape of a core can be described by the eigenval¬ 
ues of the moment of inertia tensor = f pXjXjd^x 


mass Mij = f pxipxjd^x (GOll). Let a,b and c be the 
lengths of the principal axes and a > b > c. A prolate 
core has b/a = c/a, and an oblate core b/a = 1. The 
core aspect ratio distribution in Figure [T6| shows several 
interesting features. First, most cores are prolate and 
triaxial, and only a small fraction are oblate. If cores 
form via fragmentation in isothermal filaments, they are 
expected to be initially prolate, since the fragmenta¬ 
tion scale along filaments is expect ed to have separations 
larger than the filame nt diameter (|Nagasawa|1987| |Inut- 
suka & Miyama|1992 ). In fact, the typical ti-core aspect 


ratios are smaller than the ratio of spacing to diame¬ 
ter for the fastest growing mode of filament fragmen¬ 
tation, again suggesting that cores are not exclusively 
formed as instabilities in filaments. Second, the fraction 
of cores that are prolate is higher in tcoircores than ti- 
cores, implying more evolved cores tends to be more pro¬ 
late, whereas cores formed via gravitational instability of 
filaments become more spherical in time. Last, among 
all ti-cores, more are prolate for higher Mach numbers. 


consistent with those cores in high Mach number simu¬ 
lations being slightly more evolved, as also indicated in 
the left panels of Figure [T^ 

3.3.4. Core Mass Function 

The core mass functions (CMFs) of ti-cores are shown 
in Figure [T7| To compare the CMFs from our simula¬ 
tions to the observed CMF,we match the median mass, 
and only compare the shape of different distributions 
with a logarithmic horizontal scale. This is because the 
mass scale Mg in our simulations can be adjusted by 
changing po (see Equation @)- Figure [Tt] includes the 


(shifte d) fit to full Aquila core sample oFJAndre et al. 


(20101, a nd the (shifted) fit to the gravitationally bound 


portion (Kdnyves et al. 2010). Although our statistics 
are limited for high-mass cores, for reference in Figure 
[m we also show the slope of the Salpeter mass function 
dN/dM oc 

The CMF of the full sample of ti-cores in Figure [T7| is 
similar to the shape of the observed CMF near the peak. 
Among the ti-cores, the majority of the small cores are 
weakly self-gravitating transient structures. If we only 
include the gravi tatio nally bound ti-cores (shaded his¬ 


togram in Figure 17), the distribution is simil ar to the 


(e.g., Cammie et al. 

2003 

Nakamura & Li|2008)), or the 

covariance matrix oi 

ottsets with respect to the center ot 


gravitationa lly bound cores observed in Aquila (Konyves 
et al. |2010 ). Although statistics for our siinulations 
(and observed) cores are too poor to fit a high-end mass 
slope with confidence, the distribution of high-mass ti- 
cores appear to be consistent with the log-normal fit of 
observed cores, and shows no evidence of an extended 
power-law tail. 

Figure[l8|shows the mass distribution of tcoii-cores. We 
also showTfor each model, a log-normal fit to the distri¬ 
bution. In addition. Figure includes for comparison 
the bound core Aquila fit andthe Kroupa IMF (shifted 
so that median masses match). As can be seen in Figure 
[TSl the CMFs of tcoircores show a statistically significant 
3tT in the high mass bins) deficit of high mass cores 
compared to the shape of the observed IMF, and can be 
well-fitted by relatively narrow log-normal distributions 
with standard deviations of ^ 0.3 —0.4. This is simliar to 
the bound core CMF log-normal fit in the Aquila cloud, 
which has a standard deviation 0.3. 

Figure [T^ shows a comparison between the -cores and 
the tcoii-cores, for each value of A4. It is interesting that 
the CMFs of ti-cores and Coii-cores do not match in our 
simulations. The ti-cores show a wider distribution with 
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Fig. 17.— Core mass function (CMF) of fi-cores from different Mach number simulations. To improve statistics, models of M02L4N256, 
M04L4N256, M16L1.5N512 are included in addition to the bold-face models in Tablel^ In each panel, the histogram shows the CMF of all 
fi-cores, with the shaded area denoting the gravitationally bound portion (see FigureThe black vertical line marks the resolution limit 
similar to Figure]^ The red solid and dash-dotted curves plot the fit to the observed (JiVlF (log-normal di stribution with standard deviation 
of 0.4 3 and 0.30) of the Aquila region of the full sample and the gravitationally bound portion of cores ||Andre et al.|2010| |Kdnyves et al.| 
[2 oTo| i, shifted horizontally to match the position of median core mass in the corresponding histogram, and normalized so that the total 
mass > O.OlMo (similar to their observational resolution limit) is the same as in the histogram. The median core mass for the full sample 
and the gravitationally bound portion of ti-cores are Mcoib/Mq =1.0, 0.8, 0.6, and 0.7, and Mcotb/Mq =3.7, 1.7, 1.4 and 1.1 for M =2, 
4, 8, and 16, using Mq = 23Mq in Table^ The blue dashed line segment shows the Salpeter slope dN/dM oc for reference. The 

error bars show v/iV, where N is the number of cores in each bin. 


more low-mass cores, and the characteristic core mass is 
a factor of ~ 2 smaller than that of tcoircores. However, 
the bound portion of the ti-core distribution is more sim¬ 
ilar to the tcoii-cores. 

The differences between ti-cores and tcoii-cores has in¬ 
teresting implications for interpreting observations. Ob¬ 
servations of CMFs represent a snapshot in time of 
GMCs, and often do not distinguish between gravitation¬ 
ally stable (weakly bound) or unstable (strongly bound) 


cores (e.g., Nutter fc Ward-Thompson|2007 Andre et al. 


20101. Our results suggest that it the cores dehned in 
observations are comparable to our ti-cores, then a large 
fraction of the low mass cores in observed CMFs may not 
end up collapsing, or gain a significant fraction of mass 
before they collapse. This also implies that the observed 
similarities between the CMF and IMF cannot be simply 
interpreted as a fundamental correlation between stellar 
mass and core mass, with each core contributing the same 
proportion of its mass to the final star, because we know 
that there is not a one-to-one mapping between cores ob¬ 
served at a given instant and cores at the time of collapse 
of star formation. However, the similarity between our 
bound ti-cores and our CoU-cores suggests that an analo¬ 


gous selection of just gravitationally bound cores can be 
obtained in observations, to probe the CMFs that better 
represent the initial mass reservoir for star formation. 

The difference between the CMFs of tcoii-cores in our 
models and the observed IMF could have several pos¬ 
sible reasons. First, the idealizations or physics in our 
simulations may be unfavorable for high mass core for¬ 
mation. For example, the presence of magnetic field, 
or a larger amplitude turbulent velocities might be able 
to provide additional magnetic or turbulence support, 
encouraging high mass core formatiorj^ Indeed, in our 
high-amplitude perturbation tests, 3 among the 21 cores 
are very massive with Mcoib/Mq >0.2. These high-mass 
cores are all turbulent supported with E/. /~ 2 (most 
of the smaller cores have E^/Eth < !)■ This is sugges¬ 
tive that massive turbulent supported cores may be more 
prevalent in systems with high velocity dispersion, but 
a more systematic study using global rather than local 
models is required. 

Alternatively, the final mass of a star may not be di- 

^ H owever, the magnetized simulations of |Chen &: Ostriker| 
||2014|l show no difference in median core masses from unmagne- 
tized models 
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Fig. 18.— Histogram shows core mass function (CMF) of icoii-cores, similar to Figure [T7| In each panel, the black solid curve shows the 
log-normal fit to the CMF, giving standard deviations of 0.29, 0.36, 0.39 and 0.35 for 4, 8, and 16. The median tcoll-core masses 

are Mcotb/Mq =1.9, 1.8, 1.5 and 1.3, for A4 =2, 4, 8 . and 16. usin g Mg = 23Mq in Table [l] Also shown in blue dashed curve a nd red 
dash-dotted curve, for reference, is the Kroupa IMF (|Kroupa|2001f and the bound core CMl^fit from Aquila ||K6nyves et al.|2010ll, both 
shifted horizontally to match median masses. 


rectly set by the mass of the associated tcoii-core. For 
example, more massive cores may be able to hold onto 
a larger fraction of their tcoircore mass, or capture more 
mass via competitive accretion at late times from the mu¬ 
tual reservoir. Those factors cannot be directly assessed 
in our simulations, since the sink particles are artificially 
merged, and feedback processes are not included. How¬ 
ever, there is some evidence that the sink particles in 
higher mass cores accrete faster, which is discussed in 


3.4. Sink Particles 


3.4.1. Accretion Rates of Sink Particles 


When a sink particle is created, gravity dominates over 
gas pressure in the envelope, and the gas accelerates into 
the sink, approaching free-fall. A spherical LP density 
profile with gas free-falling to the center would have a 
constant mass accretion rat43 


d{Msp/Mo) 

d{t/to) 


(27) 


In fact, pressure reduces the accretion rate somewhat 
(see below), but infall still progresses from the inside to 
the outside of the envelope. After this inside-out collapse 


® The accretion rate unit, Mg/tg = nc^/G, is equal to 5.1 X 
10“®A4Qyr“l, for Cs = 0.19km/s. 


reaches out to where the density profile flattens out and 
there is more turbulent and thermal s uppo rt, the accre¬ 
tion rate is expected to change. Figure [20| illustrates the 
accretion history of sink particles. On average, the accre¬ 
tion rate remains roughly constant, with a slight rising 
trend after the sink mass equals the tcoircore mass. How¬ 
ever, we see no clear separation between the accretion of 
the initial tidally-bound core and later stages; this is not 
surprising, as the angle-averaged density profile is con¬ 
tinuous at the effective core radius. Noticeably, there are 
outbursts in accretion rates with the peak several times 
higher than the average for a period of ~ O.Olto- This 
is due to the clumpy density structures that develop and 
fall into the sink. In reality, the infall onto the proto¬ 
star has to be processed through an accretion disk which 
is not resolved in our simulations. The additional mass 
falling into the star during outbursts of this kind may 
trigger gravitational or magnetic instabilities in the in¬ 
ner disk, leading to more dramatic outbursts in YSOs, 
for instance, the episodic accretion of F U Ori-type stars 
( Audard et al.|[20M Ohtani et al.| 2014). 

F igure |2U| also shows a trend that the average accre¬ 
tion rate increases with the co rres ponding tcoii-core mass, 
which is quantified in Figure pH The average Mgp de¬ 
pends linearly on logMcore, buT there is also a big dis¬ 
persion with a factor of ^ 2 in accretion rates at a given 
core mass. Another way to show the dependence of ac- 
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Fig. 19. — Comparison of CMFs of ii-cores and fcoii-cores (see also Figure [T7| and [T^ . In each panel, the solid histogram plots the CMF 
of tcoii-cores, with error bars showing \/]V, where N is the number of cores in each bin. The dashed histogram shows the CMF of ti-cores, 
and the shaded part the gravitationally bound portion. 


cretion rate on core mass is to look at the mass of sink 
particle at the free-fall timescale of the core. Figure 22 
shows sink particle masses at 1, 2, and 4 times the initml 
free-fall time of the tcoii-cores. The initial mass in the 
core is accreted by the sink at ^ t//, after which the 
accretion rate drops for low mass cores while it rises for 
high mass ones. At the ratio of sink mass to fcoir 
core mass is a factor of ~ 2 larger for the highest mass 
cores than the lowest mass ones. This trend provides 
one possible explanation of the discrepancy between the 
CMF in our simulations and the observed IMF, although 
in addition to this effect, the final stellar mass would de¬ 
pend on feedback and mass loss at later stages of star 
formation, which is beyond the scope of this paper. 

The sink particle accretion properties of our high- 
amplitude velocity perturbation simulations are similar 
to the low-amplitude cases discussed above. 


3.4.2. Star Formation Rate 

The total mass in sink particles, as shown in the left 
panel of Figure |23[ increases super-linearly with time, 
roughly as Mgp^totai oc This is because both the mass 
of individual sinks and the total number of sink p arti- 
cles increases linearly with time. Lee et al. (20141 also 
obtained a similar result that the accreted mass grows 
as in turbulent MHD simulations. For higher Mach 
numbers, we find that there is more mass accreted at a 
given time, both because sink particles form faster, and 
more sinks form per unit area. We note that this only 


represents the early stage of star formation when gravity 
dominates, and the surrounding gas is abundant. The 
star formation rate would have to drop at a later times 
when stellar feedback becomes important and gas in the 
CMC is exhausted or removed. 

The star formation efficiency (SFE) is shown in the 
right panel of Figure The SFE is defined as the ratio 
of total mass in all sinks to 2M.CsL^t, the total mass 
that has entered the box (or the shocked gas). As the 
mass in sinks grows oc the SFE increases linearly with 
t. Similar to all simulations without feedback (eg. iLee 
et al.| (|2014|); |Smith et al.| (|2009|)), gas is accreted on the 
free-fall timescale, resulting m a high SFE of ~ 0.3 — 0.6. 


4. SUMMARY 

In this work, we have used three-dimensional numeri¬ 
cal simulations to explore the physics of pre-stellar core 
formation and the subsequent core collapse that leads 
to the formation of protostars. We focus on local re¬ 
gions (~ Ipc) within GMCs where a supersonic turbu¬ 
lent flow converges, triggering star formation in the post¬ 
shock layer. Our simulations adopt an idealized isother¬ 
mal equation of state, include gravity from both gas and 
sink particles, and explore a wide range of Mach num¬ 
bers (A4 = 2 — 16) for the inflow gas. We terminate our 
simulations at late times, when most of the gas in the 
post-shock region has been accreted by the sink parti¬ 
cles. 

For a fixed Mach number, we have carried out a series 
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Fig. 20. — Accretion rates of sink particles for different Mach numbers. We draw a random subset of 40 sink particles in each Mach 
number case for demonstration. Since mass in the control volume around a sink is added to the sink instantaneously at the time of its 
creation, we plotted from time tcoii + i//(pLp(3Ax)), approximately the time when the mass within two times the radius of the control 
volume has fallen into the sink, until one sink particle merges with another, or the end of the simulation. The accretion rates are calculated 
using the slope of local cubic spline interpolation of the sink mass at different times, to smooth out noise from discrete time intervals (the 
time resolution of simulation output is O.OOlfo)- The color of lines show the corresponding tcow-coxe mass. The sphere on each line marks 
the time when the sink mass equals the tcoii-core mass. 
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Fig. 21 .— Average accretion rates of sink particles from + 
tff{pLp{3Ax)) to the time when the sink mass equals the cor¬ 
responding tcoii-core mass. The dotted line shows a log-linear 
fit d{Msp/Mo)/d{t/to) = 4.1 log(Mcore/Mo) + 8.5. With Cs = 
0.19km/s, The accretion rate unit is 5.1 X lO“®M 0 yr“'^, such that 
most are in the range (1 — 5) X lO“®M 0 yr“^. 


of simulations with different resolutions and box sizes to 
study potential numerical effects, especially the conver¬ 
gence of isothermal fragmentation. We have obtained the 
properties of > 200 cores in each model, for both cores 
identified at the time of individual collapse (tcoii-cores) 
and from a snapshot at the time of first singularity for¬ 
mation in each simulation (fi-cores). We study the sta¬ 
tistical properties of cores, and compare the CMFs in our 
simulations to the observed CMF and IMF. In addition, 
we trace the mass accretion rate of sink particles, and its 
dependence on the core mass. 

Our main conclusions are as follows: 


1. Cores and self-gravitating filaments form and 
evolve at the same time. Isothermal filaments do 


not collapse to become singu lar spindles (cf. Inut- 


suka & Miyama 1992, 1997) because both small- 


scale and large-scale over-densities are seeded by 
compression from turbulence, and filaments and 
cores grow simultaneously. 


2. At the time of singularity formation, all fcoii-cores 
approach the Larson-Penston asymptotic solution, 
giving a mass-radius relation of M oc r (see Fig- 
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Fig. 22. — The sink particle mass compared to the corresponding icoii-core mass at 1, 2, and 4 times the free-fall time for each core. The 
black lines plot Msp equals to 1, 2, and 4 times Mcore, and are not fits to the data. Sinks accrete close to the tcoll-core mass each free-fall 
time, although high-mass cores begin to accrete faster than low-mass cores over time (see right panel). 



Fig. 23. — Left: total mass in sink particles per unit area projected along the z direction of the simulation box. Each line represents one 
simulation run. The black solid line plots the slope of M^p total Right: the star formation efficiency (SFE, see text for definition) of 

each simulation run. 


ure 


131, with the LP density profile extending be- 


yonoseveral times the angle-averaged tidal radius 
of the cores. Most tcoii-cores are prolate, consistent 
with filament fragmentation. The mass-radius rela¬ 
tion of ti cores implies a flatter density profi le with 


& Richer 

2010 

Kirk et al. 

2013 


are only weakly bound. 


Many ti-cores 


3. At the stage of core collapse, there is a well-defined 
characteristic mass scale for fragmentation in our 
turbulent, self-gravitating simulations. The char¬ 
acteristic tcoii-core mass in our simulations con¬ 
verges as we increase the resolution and density 
threshold for sink particle creation. This mass is 
very close to the mass of a critical Bonner-Ebert 
sphere (see Equation 0)y taking the “edge” as 
equal to the mean pressure in the post-sho ck l ayer 
formed by the converging flow (See Figure 111. 


4. The characteristic core mass for the present sim¬ 
ulations is insensitive to the Mach number of the 
large-scale converging flow. As the inflow Mach 
number increases, the post-shock density and pres¬ 
sure do not vary much, due to instabilities that 


disperse the gas in the post-shock layer, leading to 
similar core properties across different Mach num¬ 
bers. For all models, the characteristic core mass is 
^ IMq. We note, however, that inclusion of mag¬ 
netic fields reduces these instabilities, such that 
the (magnetically-dominated) post-shock pressure 
is higher, and characteristic core mass is lower in 
models with higher Mach numbers (Chen & Os- 
triker 2015, in preparation!^ 


5. The CMFs of tcoii-cores in our simulations show 
a log-normal distr ibu tion with standard deviation 
0.3—0.4 (Figure 18). This is simil ar to the bound 


core distribution width for Aquila (Konyves et al. 
2010). There is a significant deficit of high mass 
cores (> 7Mq) compared to the shape of the ob¬ 
served IMF. However, the CMFs of ti-cores (Figure 


17) show a wider distribution with more low mass 


(e.g., Alves et al. 12007 

Nutter & Ward-Thompson 

2007 Andre et al. 

201 

'This suggests that ob- 


However, in the magnetized case, the characteristic core mass 
is still comparable to the critical Bonner-Ebert mass at the mean 
post-shock pressure 































































23 


served CMFs, similar to our ti-cores, may include 
cores at an early stage of evolution that will gain 
a significant amount of mass before collapsing, or 
transient structures from turbulence perturbations 
that may never collapse. Therefore, the observed 
CMF from a cloud “snapshot” may not represent 
the distribution of mass reservoirs for protostar col¬ 
lapse. We find, however, that the <i-cores that ex¬ 
ceed the critical BE mass have a distribution that 
is closer to that of the tcoircores; a similar selection 
could be made for observed cores. 


6 . The mass accretion rate of sink particles increases 
weakly with the corresponding tcoU-core mass (Fig¬ 


ure 221 , suggesting that competitive accretion may 
play a role in filling out the high-mass portion 
of the IMF. Individual sink particles accrete at 
roughly constant M^p, and continue to gain more 
mass even after they have exceeded the mass of the 
core when it first collapsed. There are outburs ts in 
« (1 — 5) X lO“®M 0 yr“^ (Figures 


20 


211 due 


to clumpy density structures forming and tailing 
into the sink. 


7. The star formation rate (measured as the total 
mass in sink particles) in each converging-flow sim¬ 
ulation increases with time as SFR oc t, as both 
the formation rate of sink particles (i.e., rate of 
cores surpassing the collapse threshold), and the 
accretion rate of individual sinks are roughly con¬ 


stant. This result is similar to the “accelerating 
star form ation” found in some other recent sim¬ 
ulations (Lee et al. 2014). For a realistic cloud, 
however, the duration of any given large-scale con¬ 
verging flow would be less than the crossing time of 
the cloud, and feedback would limit the accretion 
within any individual core. 


The present suite of models extend previous simulation 
studies of gravoturbulent fragmentation in several impor¬ 
tant ways. Most notably, by systematically varying the 
Mach number, we are able to test the influence of a key 
environmental parameter, and by systematically varying 
the resolution, we are able to demonstrate that evolution 
robustly includes a core-collapse stage with well-defined 
mass. The physical ingredients of our simulations are 
limited, however, in particular lacking magnetic fields, 
feedback from forming stars, and effects of environmen¬ 
tal evolution. Furthermore, global cloud simulations are 
required in order to incorporate the full spectrum of tur¬ 
bulence at realistic amplitudes. Other work currently 
underway will address these limitations while maintain¬ 
ing a focus on systematic surveys of the environmental 
parameter space. 
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